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Abstract  The  transport  patterns  of  non-thermal  H"*"  and 
O'*"  field-aligned  flows  from  the  dayside  cusp/cleft,  asso¬ 
ciated  with  transverse  heating  by  means  of  wave-particle 
interactions  and  in  combination  with  the  poleward  motion 
due  to  the  magnetospheric  convection  are  investigated.  This 
has  been  accomplished  by  developing  a  steady-state,  two- 
dimensional,  trajectory-based  code.  The  ion  heating  is  mod¬ 
elled  by  means  of  a  Monte  Carlo  technique,  via  the  process 
of  ion  cyclotron  resonance  (ICR),  with  the  electromagnetic 
left-hand  circular  polarized  component  of  a  broad-band,  ex¬ 
tremely  low-frequency  (BBELF)  turbulence.  The  altitude  de¬ 
pendence  of  ICR  heating  from  1000  km  to  3  Earth  radii  (Re) 
is  modelled  by  a  power  law  spectrum,  with  an  index  a,  and 
a  parameter  u^o  that  is  proportional  to  the  spectral  density 
at  a  referenced  gyrofrequency.  Because  of  the  finite  latitu¬ 
dinal  extent  of  the  cusp^left,  the  incorporation  of  the  hori¬ 
zontal  convection  drift  leads  to  a  maximum  residence  time 
to  of  the  ions  when  being  energized.  A  large  set  of  simula¬ 
tions  has  been  computed  so  as  to  study  the  transport  patterns 
of  the  H"*"  and  bulk  parameters  as  a  function  of  to,  a, 
and  1^0.  Residence  time  effects  are  significant  in  O'*"  density 
patterns  while  negligible  for  When  comparing  the  re¬ 
sults  with  analytical  one-dimensional  theories  (Chang  et  al, 
1986;  Crew  et  al.,  1990),  we  find  that  mean  ion  energies  and 
pitch  angles  at  the  poleward  edge  of  the  heating  region  are 
slightly  influenced  by  to  and  may  be  used  as  a  probe  of  ICR 
parameters  (of,  wo)-  Conversely,  poleward  of  the  heating  re¬ 
gion,  upward  velocity  and  mean  energy  dispersive  patterns 
depend  mainly  on  to  (e.g.  the  magnitude  of  the  convection 
drift)  with  latitudinal  profiles  varying  versus  to-  In  short,  the 
main  conclusion  of  the  paper  is  that  any  triplet  (to,  ct,  wq) 
leads  to  a  unique  transport-pattern  feature  of  ion  flows  asso¬ 
ciated  with  a  cusp/cleft  ionospheric  source.  In  a  companion 
paper,  by  using  high-altitude  (1.5-3  Re)  ion  observations  as 
constraints,  the  results  from  the  parametric  study  are  used  to 
determine  the  altitude  dependence  of  transverse  ion  heating 
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during  a  significant  number  of  passes  of  the  Interball-2  satel¬ 
lite. 

Key  words.  Magnetospheric  physics  (auroral  phenomena)  - 
Space  plasma  physics  (numerical  simultation  studies;  wave- 
particle  interactions) 


1  Introduction 

The  energization  and  outflow  of  ionospheric  ions  at  auro¬ 
ral  latitudes  have  been  an  area  of  active  research  over  the 
past  three  decades,  since  it  was  first  observed  by  Shelley 
et  al.  (1972).  Many  observations  and  theories  pointed  out 
that  the  ion  outflow  might  be  of  sufficient  strength  to  sup¬ 
ply  the  outer  regions  of  the  magnetosphere  (see  the  reviews 
by  Andre  and  Yau,  1997;  and  Moore  et  al.,  1999).  The  en¬ 
larged  cusp/cleft  region  has  been  identified  as  a  major  source 
of  ionospheric  ions  for  the  magnetosphere  (Lockwood  et  al, 
1985a;  Thelin  et  al.,  1990).  This  region  is  located  in  the 
dayside  auroral  zone,  extending  between  about  9  and  15  h  in 
magnetic  local  time  (MLT)  with  a  latitudinal  width  of  a  few 
degrees. 

In  the  cusp/cleft,  much  of  the  ionospheric  ion  outflow  is 
directly  caused  by  energization  of  the  major  ion  species  (H"^, 
O'*')  transverse  to  the  geomagnetic  field  (Andre  and  Yau, 
1997,  0ieroset  et  al.,  2000).  In  these  regions,  the  heated 
distributions  form  so-called  conics  in  velocity  space,  with  a 
peak  flux  in  an  oblique  direction  to  the  geomagnetic  field. 
Statistical  studies  of  the  evolution  of  ion  conics  versus  alti¬ 
tude  in  the  cusp/cleft,  from  1  to  4  Earth  radii  (Re),  using 
DE-1  (Peterson  et  al.,  1992)  and  up  to  10  000  km,  using  Ake- 
bono  data  (Miyake  et  al.,  1993),  pointed  out  that  their  char¬ 
acteristic  energies  increase  while  their  apex  angle  (the  cone 
angle  with  respect  to  the  upward  direction)  decreases  much 
more  slowly  than  expected  from  adiabatic  motion.  Both  of 
these  findings  are  consistent  with  gradual  heating  of  ion  con¬ 
ics  as  they  move  upward  along  the  geomagnetic  field  lines. 
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Depending  on  the  altitude  of  observations,  typical  tempera¬ 
tures  of  ion  conics  range  from  10  eV  up  to  at  least  a  few  keV 
(Moore  et  al,  1999). 

In  the  dayside  polar  cap  poleward  of  the  cusp/cleft,  obser¬ 
vations  by  DE-1  revealed  that  outflowing  ions,  first  heated 
in  the  cusp/cleft,  overcome  gravity  via  the  mirror  force  and 
drift  poleward  under  the  effect  of  the  magnetospheric  con¬ 
vection  to  form  the  so-called  ion  fountain  (Lockwood  et  al, 
1985b).  A  scheme  of  a  satellite  moving  across  the  ion  foun¬ 
tain  is  given  in  Fig.  1.  When  observed  at  high  altitudes, 
outflowing  ion  distributions  are  conical  at  high  energies  in 
the  heating  region,  and  field-aligned  at  lower  energies  pole- 
ward  to  the  heating  region  (Horwitz,  1986;  Knudsen  et  al, 
1994;  Dubouloz  et  al,  1998).  Qualitatively,  this  picture  can 
be  formed  by  ion  transverse  heating  in  a  region  of  finite 
latitudinal  extent,  followed  by  adiabatic  convective  flow  to 
the  satellite  orbit.  Poleward  of  the  heating  region,  the  lat¬ 
ter  stage  contributes  to  a  velocity-filter  effect  as  evidenced 
by  the  latitude-energy  dispersion  on  ion  data.  Many  models 
based  on  ion  trajectory  calculations  have  been  developed  in 
order  to  demonstrate  that  the  high-altitude  observed  outflows 
originate  from  the  cusp/cleft  (Horwitz  and  Lockwood,  1985; 
Knudsen  et  al,  1994;  Dubouloz  et  al,  1998,  2001).  How¬ 
ever,  the  effect  of  non-adiabatic  processes,  such  as  transverse 
heating,  have  not  been  incorporated  when  modeling  the  ion 
transport  in  the  cusp/cleft. 

Because  all  ionospheric  species  achieved  non-thermal  en¬ 
ergies  when  observed  in  the  cusp/cleft,  the  heating  is  caused 
by  the  essentially  perpendicular  component  of  electric  fields 
of  oscillating  within  some  frequency  range  (Andre  and 
Yau,  1997).  Recent  statistical  observations  from  the  Freja 
(Norqvist  et  al,  1998)  and  FAST  (Lund  et  al,  2000)  satel¬ 
lites  pointed  out  that  a  major  part  of  ion  transverse  heat¬ 
ing  events  are  associated  with  an  enhancement  of  the  broad¬ 
band,  extremely  low-frequency  (BBELF)  wave  turbulence. 
The  BBELF  turbulence  covers  frequencies  from  less  than 
1  Hz  up  to  several  hundred  Hz,  thus  including  the  H*^  and 
gyrofrequencies  at  altitudes  from  about  1000  km  up  to 
a  few  Rf.  Its  profile  versus  frequency  exhibits  generally  a 
power  law  spectrum  (Kintner,  1976;  Gumett  et  al,  1984), 
which  may  be  modelled  by  an  index  a.  In  terms  of  ion  en¬ 
ergization,  the  theory  of  ion  cyclotron  resonance  (ICR)  heat¬ 
ing  by  the  electromagnetic  left-hand  circular  polarized  com¬ 
ponent  around  the  ion  gyrofrequencies  (Chang  et  al,  1986) 
seems  the  most  plausible  mechanism.  Adopting  a  power  law 
altitude-dependent  ion  heating  profile  and  computing  Monte- 
Carlo  simulations,  many  authors  confirmed  that  only  a  small 
fraction  of  the  wave  intensity  in  the  BBELF  spectrum  is 
needed  to  cause  the  observed  ion  energies  through  ICR  heat¬ 
ing  (Retterer  et  al,  1987;  Crew  et  al,  1990;  Andr^  et  al, 
1990;  Norqvist  et  al,  1996).  Except  for  the  work  of  Andre 
et  al.  (1990),  where  detailed  wave  and  particle  observations 
available  for  one  event  have  been  used  to  carry  out  a  2-D 
Monte  Carlo  simulation,  all  the  previous  models  neglected 
the  horizontal  convection  drift.  Because  of  the  finite  latitudi¬ 
nal  extent  of  the  cusp/cleft,  the  effect  of  a  poleward  convec¬ 
tion  leads  to  a  limited  residence  time  tp  of  ions  when  being 


Fig.  1.  Sketch  illustrating  the  transport  patterns  of  ion  outflows  from 
the  cusp/cleft  in  conjunction  with  a  satellite  crossing.  The  trajecto¬ 
ries  are  projected  in  a  meridian  plane  and  are  functions  of  invariant 
latitude  and  altitude. 


heated. 

In  the  present  study,  a  steady-state,  2-D,  Monte  Carlo 
model  has  been  developed.  A  large  set  of  simulations  is  com¬ 
puted  so  as  to  study  the  effects  of  different  geophysical  pa¬ 
rameters,  such  as  the  altitude  dependence  of  the  ion  heating 
rate  or  the  magnetospheric  convection  drift,  on  the  transport 
patterns  of  H"*"  and  O’’"  ions.  The  aim  of  this  parametric  study 
is  focussed  on  determining  the  altitude  dependence  of  ion 
transverse  heating  for  a  specific  event,  providing  that  high- 
altitude  (below  3  Re)  ion  observations  inside  and  poleward 
of  the  cusp/cleft  are  available. 

The  outline  of  the  paper  is  as  follows.  A  description  of  the 
simulation  model  is  given  in  Sect.  2,  along  with  a  summary 
of  the  key  parameters  used  as  input  to  a  simulation.  For  a 
few  simulations,  brief  results  and  some  tests  are  discussed  in 
Sect.  3,  The  results  of  the  parametric  study  applied  to  H+ 
and  ion  flows  are  discussed  quantitatively  in  Sect.  4.  A 
summary  of  conclusions  is  presented  in  Sect.  5. 
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Fig.  2.  Transport  patterns  of  0+  moments  in  the  (X,  s)  plane  :  (a)  density  in  cm~^  for  a  value  of  5  000  cm”^  at  the  low-altitude  boundary, 
(b)  upward  mean  velocity  in  km  s”  ^  (c)  mean  energy  in  eV,  (d)  fluence  in  Log(cm"2  s"  ^ ),  (e)  perpendicular  mean  energy  in  eV,  (f)  parallel 
mean  energy  in  eV,  (g)  mean  pitch  angle  in  and  (h)  energy  density  flux  in  Log(/>tWm"2).  The  vertical  white  line  corresponds  to  the 
poleward  heating  boundary  (X  =1). 
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2  Model 

2. 1  Model  geometry  and  applicability 

The  model  geometry  considered  here  is  similar  to  the  one 
used  by  Andre  et  al.  (1990),  According  to  the  scheme  in 
Fig.  1,  the  coordinate  system  is  defined  by  two  spatial  di¬ 
mensions:  s  and  X,  where  s  denotes  the  length  of  the  ge¬ 
omagnetic  field  line  from  its  footprint,  plus  the  Earth’s  ra¬ 
dius.  For  high-latitude  regions  associated  with  opened  field¬ 
lines  and  in  our  altitude  range  (.y  <4  Re),  the  curvature  of  the 
magnetic  field  lines  may  be  neglected  and  s  is  close  to  the 
geocentric  distance.  On  the  other  hand,  we  take  into  con¬ 
sideration  the  divergence  of  the  magnetic  field  lines  by  in¬ 
troducing  an  altitude  dependence  proportional  to  l/.y^.  The 
dayside  heating  regions,  such  as  the  cusp  or  the  cleft,  have 
longitudinal  boundaries  at  equal  latitudes  (Newell  and  Meng, 
1992).  Hence,  by  applying  a  2-D  model  to  a  data  set,  it  is 
assumed  that  the  heating  region  is  longitudinally  homoge¬ 
neous.  The  poleward  drift  Vp  caused  by  the  x  B  convec¬ 
tion  is  taken  into  consideration  by  introducing  a  poleward 
abscissa  jc,  where  x  is  a  measurement  of  the  poleward  dis¬ 
tance  from  the  equatorward  border  of  the  heating  region.  For 
convenience,  we  introduce  a  distance  X  normalized  with  re¬ 
spect  to  the  thickness  of  the  region.  The  relation  between  X 
and  the  invariant  latitude  A  is  given  by: 


X=:(A- A,^)/A,  (1) 

where  Agq  is  the  invariant  latitude  of  the  equatorial  boundary 
and  A  the  latitudinal  width  of  the  heating  region. 

In  the  present  study,  in  order  to  have  a  reasonable  number 
of  free  parameters  in  the  simulation,  we  consider  a  constant 
drift  Vp  along  the  X  axis.  In  the  (X,  s)  frame,  the  effect  of 
the  drift  leads  to  a  limited  residence  time  for  ions  given  by 
to  =  A/vp,  where  A  is  the  width  of  the  heating  region  in 
kilometers.  In  the  high-latitude  ionosphere,  Vp  ranges  from 
100  to  1000  ms~\  and  A  is  about  a  hundred  kilometers. 
Therefore,  to  varies  typically  between  100  and  1000  s. 

Ionospheric  convection  patterns  depend  strongly  on  the  in¬ 
terplanetary  magnetic  field  (IMF)  orientation  (Ruohoniemi 
and  Greenwald,  1996).  In  the  same  way,  the  position  of 
the  cusp  may  be  strongly  modulated  by  the  IMF  orientation 
(Chen  et  al.,  1997).  Hence,  when  assuming  a  constant  drift  in 
the  model,  it  implies  steady  IMF  conditions  during  the  trans¬ 
port  of  the  ion  field-aligned  flows  between  the  ionosphere 
and  the  satellite  orbit. 

The  simulation  region  extends  in  s  from  the  topside  iono¬ 
sphere  (at  1000  km  altitude)  to  4  Re,  and  in  latitude  from 
X  =  0  to  X  =  10.  As  onset  to  a  simulation,  H"*"  and  O"*"  test 
particles  are  launched  from  the  topside  auroral  ionosphere 
(s  =  .so  and  0<X<1)  and  possibly  along  the  equatorward 
edge  of  the  heating  region,  ^en  considering  homogeneous 
heating  region  in  latitude  and  a  constant  drift,  the  ions  as¬ 
sociated  with  those  two  flux  tubes  do  not  mix.  Hence,  the 
transport  of  those  populations  may  be  computed  as  indepen¬ 
dent  elements,  with  the  equatorward  population  leading  to 


a  more  accurate  description  of  the  low-latitude  side  of  the 
heating  region. 

The  structure  of  the  steady-state  ion  flow  patterns  is  carried 
out  from  the  computed  ion  trajectories.  The  ions  dynamics 
include  the  effect  of  the  wave-particle  interactions  (WPI),  the 
different  macroscopic  forces  (gravity  g,  mirror  force),  and 
the  poleward  drift  motion  caused  by  the  convection.  On  the 
other  hand,  the  model  neglects  the  curvature  of  the  field  lines 
and  therefore  their  effects:  curvature  and  gradient  drifts,  and 
centrifugal  acceleration  (Cladis,  1986).  The  effects  of  WPI 
are  modelled  using  a  Monte  Carlo  technique,  as  described  in 
the  next  section. 


2.2  Ion  heating  model 
2.2. 1  Kinetic  formulation 


For  a  type  of  ion  /,  the  time  evolution  of  its  distribution  func¬ 
tion  F/(r,  V,  0  in  the  phase  space  (r,  V)  is  given  by  the 
kinetic  equation  (Ichimaru,  1973): 


g  +  -^(£  +  VxB)l-^) 
m/  J  oV  J 


Fiir,  V)=-^  Di.  V),  (2) 

where  qt  and  rrii  are  the  ion  charge  and  mass,  respectively. 
The  right-hand  side  term  of  Eq.  (2)  expresses  the  heating 
of  ions  resulting  from  the  wave-particle  interactions,  as  de¬ 
scribed  by  a  quasi-linear  velocity  difiusion  tensor  Vi  Di .  Be¬ 
cause  Fi(r,  V,  0  is  gyrotropic  in  our  problem,  and  its  mea¬ 
surement  takes  place  on  times  greater  than  the  ion  cyclotron 
period,  we  may  consider  the  evolution  of  the  distribution 
which  is  the  gyrophase  average  of  F|(r,  V).  The  phase 
space  velocity  is  then  decomposed  in  two  components  V\\ 
and  Vx,  parallel  and  perpendicular  to  B,  respectively.  Con¬ 
sidering  the  approximations  discussed  above  (i.e.  neglecting 
the  field-line  curvature,  including  the  convection  drift  only), 
the  evolution  of  fi(s,  X,  V(|,  Vx)  is  given  by  the  equation 
L Ms,  X,  Vi|,  Vx)  =  0  (Crew  and  Chang,  1985),  where  L 
is  a  kinetic  operator,  that  includes  the  phase  space  flow,  and 
can  be  expressed  as: 

^  ds  r/)  9X  nti  ds  9V||  2  ds 


The  spatial  variations  of  ft  are  entirely  controlled  by  the  first 
two  terms  in  the  right-hand  side.  The  first  term  describes  the 
motion  along  the  magnetic  field  lines,  while  the  second  term 
((l/ri))9/9X  =  Vpd/dx)  follows  the  poleward  motion  due 
to  the  B  X  B  convection.  The  next  terms  force  the  distribu¬ 
tion  to  respond  to  the  action  of  macroscopic  forces,  such  as 
the  potential  gradient  d^/ds  =  — (F||  H-  {rnilqi)g{s))  due  to 
the  gravity  g(.y)  and  a  possible  parallel  electric  field  E\\  (third 
term)  and  the  mirror  conversion  (fourth  term).  We  remind  the 
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Fig.  3.  Transport  patterns  of  moments  in  the  (X,  s)  plane.  The  format  is  the  same  as  in  Fig.  2. 
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AiUtude  (km) 

Fig.  4.  Evolution  of  the  perpendicular  (red  points)  and  parallel  (blue 
points)  energies  of  (top)  and  (bottom)  test  particles  inside 
the  heating  region  as  a  function  of  altitude.  The  corresponding  solid 
curves  are  inferred  from  the  mean  ion  theory  of  Chang  et  al.  (1986). 

reader  that  the  mirror  term  has  this  simple,  one-dimensional 
form  because  of  the  neglect  of  curvature  and  gradient  drifts 
and  centrifugal  acceleration,  as  discussed  in  Sect.  1 .  In  the 
case  of  ICR  heating  at  the  ion  gyrofrequency  fd,  the  compo¬ 
nent  of  the  diffusion  tensor  perpendicular  to  B  is  dominant, 
i.e.  D|X  >  For  wavelengths  greater  than  the  ion  gy- 

roradius,  D/x  can  be  expressed  as  (Chang  et  al.,  1986;  Crew 
etal.,  1990): 

A.xO)  «  (4) 

where  Silfd]  is  the  left-hand  component  of  the  electric  field 
spectral  density  at  fd.  This  results  in  a  useful  simplification 
of  having  a  velocity-independent  diffusion  term. 

2.2.2  A  Monte  Carlo  technique 

The  solution  fi(s,  X,  Vn,  V±)  of  Eq.  (3)  is  carried  out  using  a 
Monte  Carlo  technique,  as  usually  applied  to  investigate  this 
kind  of  problem  (Retterer  et  al.,  1983,  1989).  In  the  high- 
altitude  cusp/cleft  region,  the  ionospheric  population  is  a  mi¬ 
nority  constituent  of  the  plasma  (when  compared  to  the  mag¬ 
netosheath  ion  population).  Therefore,  we  may  treat  them 
as  test  particles  in  externally  imposed  fields,  rather  than  cal¬ 
culating  the  fields  self-consistently.  From  the  initial  source 


Fig.  5.  Evolution  of  the  integrated  particle  (left-hand  side)  and  the 
energy  density  (right-hand  side)  fluxes  of  O'*"  (top)  and  (bot¬ 
tom)  ions.  Red  and  blue  curves  correspond  to  the  integrated  fluxes 
over  the  poleward  border  (X  =  1)  up  to  an  altitude  z,  and  over 
the  horizontal  extent  of  the  heating  region  (0  <  X  <  1)  at  an  al¬ 
titude  z,  respectively.  The  black  curves  correspond  to  the  sum  of 
all  contributions  and  the  dashed  curve  (right-hand  side)  denotes  the 
transferred  from  the  waves  to  particles  below  z.  All  quantities  are 
normalized  to  the  values  at  the  low-altitude  boundary. 


distribution  5'o(r,  V)  in  velocity  and  space,  the  calculation 
of  fi(Sy  X,  V||,  V±)  proceeds  by  following  the  kinetic  gyro- 
centers  of  test  particles.  The  motion  of  an  ion  gyrocenter  is 
described  by  the  set  of  equations: 

X  =  \/tD  (5) 

rtiiVw  =  qiE\\  -  niigis)  -  Wi(dln 5/dy),  (6) 

where  V\\  and  V±  are  the  parallel  and  perpendicular  veloc¬ 
ity  components,  and  W±  is  the  perpendicular  energy.  Equa¬ 
tion  (5)  describes  the  poleward  drift  motion  along  X  due  to 
the  convection.  Equation  (6)  gives  the  variation  of  the  par¬ 
allel  velocity  under  the  effect  of  macroscopic  forces,  as  de¬ 
scribed  in  Sect,  2.2.1.  In  the  simulations,  we  set  E\\  =  0  for 
several  reasons.  First,  introducing  E\\  ^  0  would  be  phys¬ 
ically  meaningful  only  if  its  two-dimensional  variations  can 
be  taken  into  consideration  and  in  a  self-consistent  manner 
that  imposes  the  inclusion  of  all  ion  and  electron  populations 
in  the  model  (see,  for  example,  Jasperse,  1998).  Because 
we  focus  on  a  parametric  study  associated  with  a  reasonable 


M.  Bouhram  et  al.:  Modeling  transverse  heating  and  outflow  of  ionospheric  ions 


1759 


(a) 


(b) 


(c) 


1000  10000  100000 
Log(Altitude) 


Loo(Altitude) 


Fig.  6.  Evolution  of  O'*"  moments  along  the  poleward  heating 
boundary:  (a)  density,  (b)  upward  mean  velocity,  (c)  perpendic¬ 
ular  (red)  and  parallel  (blue)  mean  energies.  The  solid  and  dashed 
curves  are  associated  with  a  source  distribution  without  temperature 
and  with  a  parallel  temperature  of  0.3  eV,  respectively. 


number  of  free  parameters  in  the  simulations,  the  introduc¬ 
tion  of  an  E\\  ^  0  component  is  beyond  the  scope  of  the 
paper.  However,  it  should  be  noted  that  in  other  regions  than 
the  dayside  cusp/cleft,  such  as  in  the  downward  current  pre¬ 
midnight  aurora,  the  field-aligned  currents  j\\  and  E\\  are  usu¬ 
ally  stronger.  Then  E\\  has  a  direct  consequence  in  the  ion  en¬ 
ergization,  since  the  heating  rate  is  found  to  be  proportional 
to  the  dissipation  term  j\\E\\  (Lynch  et  al.,  2002).  Hence,  in 
the  companion  paper,  simulations  are  consequently  applied 
to  data  sets  in  the  cusp/cleft,  where  the  effect  of  E\\  is  esti¬ 
mated  to  be  negligible,  using  the  observed  electron  and  ion 
distributions  at  high  altitude. 

In  the  plane  perpendicular  to  B,  the  effect  of  a  stochastic 
ion  heating  results  in  disturbing  the  velocity  V _i_  with  a  ran¬ 
dom  impulse  A  V  at  each  time  step  At.  The  two  components 
A  Vjc .  A  Vy  of  A  Vj_  are  generated  by  a  centered  Gaussian  dis¬ 
tribution,  such  that  (A  V^)  =  [AVy)  =  2D±At  (Retterer  et 
al.,  1983,  1989).  The  angle  0  That  describes  the  orientation 


of  AVI  with  respect  to  the  ion  gyrophase  is  assumed  to  be 
random.  Thus,  the  variation  A  of  the  perpendicular 

energy  due  to  an  impulse  in  a  time  step  At  is  given  by: 

AlVix.rfi  =  j»w/(Vx  +  AVx)^  -  ^/WilVxP 

1  9 

=  m/VxAVxcose  +  -/n,|AVxp.  (7) 

Because  6  is  randomly  oriented,  only  the  last  term  of  the 
right-hand  side  of  Eq.  (7)  contributes  statistically  to  a  mean 
ion  heating  rate  given  by: 

W/x(^)  =  2/n/Ai.  (8) 

Inside  the  heating  region  (X  <  1),  the  time  step  At  must  sat¬ 
isfy  the  condition  lV,  i(.s)  x  At  «:  W,  x(*s),  so  that  the  test 
particle  keeps  its  energy  Wi  and  first  adiabatic  invariant  /x,* 
constant  between  two  random  impulses.  Typically,  we  use 
At  =  with  s  ~  0.01.  Outside  the  heating  re¬ 

gion,  the  transverse  heating  is  turned  off,  and  At  is  imposed 
by  the  adiabatic  motion  (typically  At  0.3/0)^/)). 

2.2.3  Altitude  dependence  of  ion  heating 

As  discussed  in  Sect.  1,  the  wave  power  spectra  exhibit  gen¬ 
erally  a  power  law  in  frequency  (Gumett  et  al.,  1984;  Andre 
et  al,  1990):  S{f)  =  5o(/o//)“,  where  So  denotes  the  spec¬ 
tral  density  at  a  given  frequency  /o,  and  a  the  spectral  index. 
As  assumed  in  previous  models  (Andre  et  al,  1990;  Crew  et 

al,  1990),  the  wave  spectrum  does  not  vary  versus  altitude, 
and  the  fraction  due  to  the  left-hand  polarized  component 
is  constant  versus  frequency.  Recently,  a  statistical  profile 
has  been  reported  for  altitudes  up  to  10000  km  using  wave 
measurements  from  Akebono  correlated  with  transversely 
heated  ions  (Kasahara  et  al,  2001).  Plots  on  the  altitude  dis¬ 
tribution  of  the  wave  power  at  a  frequency  of  5  Hz,  which 
corresponds  roughly  to  the  O'*"  gyrofrequency  at  8000  km, 
pointed  out  that  the  wave  power  remains  approximately  con¬ 
stant  versus  altitude  in  the  dayside  cusp/cleft  (Plate  5b  of 
Kasahara  et  al,  2001),  ranging  from  about  9  x  10“^  to 
2.5  X  10“^  V^m“‘^Hz"^  This  latter  result  is  at  least  com¬ 
patible  with  the  assumption  that  the  wave  power  spectrum  is 
not  altitude  dependent.  Considering  these  assumptions,  since 
the  ion  gyrofrequency  is  proportional  to  1/s^,  the  altitude  de¬ 
pendence  of  the  ion  heating  rate  is  given  by: 

w, x(^)  =  X  (s/so)^",  (9) 

where  W/xf^o)  =  (9f/2»J,)Sx[(/c,(Jo))]  denotes  the  heat- 
ing  rate  at  the  lower  boundary  distance  ^o-  Por  heavier  ions 
of  mass  number  M,*  =  mi/mn,  lV/x(.so)  scales  as 
For  comparing  simulation  results  for  different  ions,  it  is  more 
convenient  to  associate  the  altitude  profile,  given  by  Eq.  (9), 
with  two  parameters  that  do  not  depend  on  the  type  of  ion.  In 
doing  so,  we  introduce  the  second  parameter  wo  given  by: 

Wo  =  '  = 
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where  u;o  has  the  dimension  of  energy.  On  the 
basis  of  typical  wave  observations  (see,  for  exam¬ 
ple,  Andre  et  al.  (1990),  where  0.5<or<2.5  and 
10“^  <  S(fci)  <  10”^  V^m"^Hz~*  for  0+  at  j  =  4  /?£,  wq 
ranges  from  0.1  to  lOOOeV.  Let  us  note  that  using  a  power 
law  dependence  versus  s  is  also  consistent  with  the  statistical 
results  on  the  evolution  of  ion  conics  along  the  field  lines, 
which  suggest  a  gradual  heating  of  ions  in  the  cusp/cleft  (Pe¬ 
terson  et  al.,  1992;  Miyake  et  al.,  1993). 

2.3  Steady-state  solution 

There  is  a  fundamental  identity  between  the  time  integrated 
solution  of  an  initial  value  problem  and  a  steady-state  prob¬ 
lem  (Tajima,  1989).  Let  us  note  F(r,  V,  t)  for  the  distribu¬ 
tion  fiinction  associated  with  test  particles  at  a  time  t  during 
the  simulation.  The  evolution  of  F  is  governed  by  the  dy¬ 
namical  kinetic  equation: 

^  +  LF  =  0.  (11) 

at 

where  L  is  the  kinetic  operator,  as  described  in  Eq.  (3).  At 
/  =  0,  the  distribution  fonction  is  determined  by  flie  initial 
condition:  F  =  5'o(r,  V)  =  niiso)8{s  —  so)fi{V)  (injection 
Bis  =  Jo)  and  after  a  time  AT,  F  =  0  (i.e.  no  more  test 
particles  inside  the  simulation  domain).  Integrating  Eq.  (1 1) 
from  f  =  0  to  r  =  AT,  with  fi  =  AT^^  Fdt,  yields: 

Ly;=o,  (12) 

with  the  boundary  condition  fi  ^  5o(r,  V)  aXs  =  jq*  Note 
that  the  solution  f  (j,  X,  V\\ ,  VI)  of  this  time-averaged  equa¬ 
tion  is  the  same  as  in  the  steady-state  kinetic  equation  for  a 
source  5o  (Eq.  3).  Therefore,  we  may  carry  out  the  solu¬ 
tion  fi(s,  X,  Vjj,  Vx)  by  integrating  over  time  the  motion  of 
test  particles  associated  with  the  source  Sq.  In  this  approach, 
each  test  particle  contributes  to  the  phase  space  density  of  a 
cell  by  the  time  spent  in  that  cell.  In  the  same  way,  the  steady 
state  moments  of  fi(s,  X,  V||,  Vx)  are  calculated.  As  output 
of  a  simulation,  only  the  following  moments  are  conserved 
at  any  location  (s,  X): 


y»V|[=:+O0  ^Vx=+CX) 

mis,  X)=  I  I  Ms,  X,  V||,  Vx)x 


7v||=-oo  7vj.=o 

2;rVxdVjdVx 

(13) 

/•V||=+oo 

«,(^.X)  =  (l/n,)  / 

7  V||  =  -00 

^Vx=4*oo 

/  (V||y;(j,  X.  V|1,  Vx)  X  2n-VxdV||dVx 

Jvj.=0 

(14) 

Ms,  X)  =  mis,  X)  X  mis,  X)  x  Ais) 

(15) 

/•’'ll  =+00 

X)  =  imi/ltii)  1 

7V||=— 00 

f^Vi^+oo 


f  ^  V(Ms.  X,  V,| ,  Vx)  X  27r  VxdV|,dVx 
JVx=0 


(16) 


^V|l=+oo 

/V|l=-oo 


Wi±(s,  X)  =  (nti/lm)  f 

Jvn 

/  yly;(i,X,V||,Vx)x2;rVxdV||dVx  (17) 

Jv±=o 

,V||=+oo  ^Vj.=+oo 

Jwi(s,X)  =  {mi/2)  /  V||(V|f+Vi) 

JV||=-oo  JV±^0 


Ms.  X,  V||,  Vx)  X  2;rVxdV|idVx.  (18) 

In  Eqs.  (13)-(18),  i  denotes  the  type  of  ion  (H"*"  or  O"^)  and 
tti,  Ui,  Jf,  Wii\,  Wi±,  and  Jwi  are  the  ion  density,  upward 
drift  velocity,  ion  fluence,  mean  parallel  energy,  mean  per¬ 
pendicular  energy,  and  energy  density  flux,  respectively.  In 
Eq.  (15),  the  parameter  A(s)  =  B(so)/B(s)  —  (s/sq)^  is  a 
geometric  factor,  due  to  the  divergence  of  the  geomagnetic 
field  lines.  If  the  convection  drift  is  turned  off,  /4(s)  would 
correspond  to  the  ratio  of  cross-section  areas  of  a  flux  tube 
between  .so  and  s.  Ait  =  0,  about  10^  and  H"*"  test  par¬ 
ticles  are  injected  with  an  upward  mean  velocity  «i||(5o)  of 
0.4  km  s”  1  and  5  km  s“^  respectively.  Such  values  are  con¬ 
sistent  with  typical  observations  at  this  altitude  (Chandler  et 
al.,  1991).  The  parallel  temperature  7}||(^o)  associated  with 
those  populations  is  about  a  fraction  of  1  eV.  In  most  of  the 
simulations,  we  set  Ti\\  (.so)  =  0.  However,  we  show  in  Sect.  3 
that  the  results  at  higher  altitudes  are  insensitive  to  the  mo¬ 
ments  characteristics  of  the  initial  distribution  5o(r,  V). 


3  First  results,  model  tests 
3 . 1  Simulation  example 

We  ran  a  Monte  Carlo  simulation  with  the  following  param¬ 
eters:  time  of  drift  tp  =  500  s,  spectral  index  a  =  \,  and 
wq  =  62  eV  (see  Eq,  1 1).  Figures  2  and  3  show  the  and 
H"’"  moments,  respectively,  in  a  plane  X  —  5.  Inside  the  heat¬ 
ing  region,  the  density  decreases  versus  s  faster  than  \/s^,  in 
order  to  satisfy  the  flux  conservation.  Conversely,  the  den¬ 
sity  decreases  as  a  function  of  X  slowly  outside  the  heating 
region  (X  >  1).  In  the  panels  showing  the  upward  mean  ve¬ 
locity  Ui  and  the  parallel  mean  energy  Wi  ||,  we  can  identify 
at  larger  altitudes  a  dispersion  versus  X  due  to  the  convec¬ 
tion.  The  spread  of  ion  flux  tubes  also  appears  clearly  in 
the  contours  of  ion  fluence  7,  .  This  effect  is  even  more  sig¬ 
nificant  for  heavier  ions,  such  as  O"^,  which  have  a  longer 
time-of-flight  for  reaching  a  geocentric  distance  s.  This  may 
also  explain  the  absence  of  H"*"  at  large  poleward  distances 
X. 

Inside  the  heating  region,  the  perpendicular  temperature 
Tx  and  the  average  pitch  angle  of  ion  conics  increase  grad¬ 
ually  versus  s.  This  effect  is  a  consequence  of  the  heat¬ 
ing  profile  increasing  versus  s,  with  the  perpendicular  ener¬ 
gies  being  progressively  transferred  into  parallel  energies  via 
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Fig.  7.  Evolution  of  O’*"  moments  along  the  poleward  heating 
boundary:  (a)  density,  (b)  upward  mean  velocity,  (c)  perpendic¬ 
ular  (red)  and  parallel  (blue)  mean  energies.  The  solid  and  dashed 
curves  are  associated  with  a  source  distribution  injected  at  an  alti¬ 
tude  of  1000  and  2000  km,  respectively. 


the  mirror  force.  Poleward  of  the  heating  region  (X  >  1), 
only  the  mirror  force  is  acting  along  with  gravity  for  heavy 
ions,  and  leads  to  a  steep  decrease  of  T±_  versus  X.  Another 
counterintuitive  result  appears  in  the  mean  pitch  angle  panels 
(Figs.  2g  and  3g),  where  and  contours  exhibit  differ¬ 
ent  values  outside  the  heating  region  at  the  lower  boundary. 
Here,  we  should  note  that  the  lower  limit  for  0“^  is  0°  (down¬ 
ward  flow)  instead  of  90°.  Unfortunately,  this  does  not  ap¬ 
pear  in  Fig.  2g,  because  the  color  scale  is  set  between  90° 
and  180°,  in  order  to  keep  the  discontinuity  at  the  poleward 
heating  boundary  (X  =  1)  visible  for  both  ions.  We  will 
see  in  Sect.  4.1.2  that  the  difference  can  be  explained  by  the 
influence  of  the  gravity  for  low-energy  heavy  ions,  with  the 
corresponding  trajectories  being  parabolic  in  the  low-altitude 


polar  cap  (Horwitz,  1984). 
3.2  Tests 


3.2. 1  Validation  of  ion  trajectories 

In  the  case  tj)  ^  oo  (1-D  model,  i.e.  no  convection  drift) 
and  neglecting  gravity,  it  is  possible  to  derive  an  asymptotic 
solution  of  the  equations  describing  the  motion  of  a  mean 
particle  (Chang  et  al.,  1986).  In  the  limit  of  high  energies, 
the  perpendicular  and  parallel  energy  components  (W±,  Wy) 
are  obtained  by  using  the  following  expressions: 

J  X  2a+2/3 

*'-W=2./>[(3»+i;(6.  +  ll)F»U 


,  .  ^  (6o!  +  2)u;oMf°  /^x  2«+2/3 

^  2V3[(3Q;  +  l)(6a-hll)]2/3V^0>' 


(20) 


Figure  4  shows  the  evolution  of  Wi±  and  Win  fimctions 
of  s  for  wq  =  62eVs”*  and  a  =  1.  The  convection  drift 
has  been  turned  off  so  that  the  test  particle  stays  inside  the 
heating  region.  Here,  we  found  that  the  evolution  of  W/j_ 
and  Will  is  in  good  agreement  with  the  mean  particle  theory 
developed  by  Chang  et  al.  (1986).  These  results  allow  us  to 
validate  the  computation  of  ion  trajectories  inside  the  heating 
region. 


3.2.2  Global  conservation  laws 


For  any  simulation  we  run,  it  is  possible  to  validate  the  spatial 
structure  of  ion  flows  by  using  global  conservation  laws.  By 
integrating  the  kinetic  equation  L//  =  0  over  all  velocity 
space,  we  obtain  the  local  conservation  equation  of  particle 
flux  densities: 


1  BNj 
to  dx 


=  0, 


(21) 


where  Ni  =  A(s)n/(.y,  X).  When  integrating  Eq.  (21)  over 
the  box  limited  by  0  <  X'  <  X  and  5o  <  we  obtain 

the  global  conservation  law  of  particle  flux  densities: 


r  '  Ni(X,s')to^ds'+  r  ^  Ni(X',s)ui(X',s)dX' = 
Js’=so  Jx'=0 

r  Ni(0,  r  Ni(X',  so)uiiX',  io)dX'.(22) 

Js'=So  JX'=r0 


By  multiplying  the  kinetic  equation  by  W  =  w,-  ( + V|j^)/2, 
and  by  integrating  over  all  velocity  space,  we  obtain  the  local 
conservation  equation  of  energy  flux  densities: 

tD  dX 


i-(Ni(Qi  +  Uiqi<Pi))  =  N,  WiAs),  (23) 

where  qi<t>i  includes  the  gravitational  and  electric  potential 
energy,  and  Qi  denotes  the  energy  flux.  Equation  (23)  may 
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Fig.  8.  Evolution  of  0“^  (left)  and  H"*"  (right)  moments  along  the  poleward  heating  boundary.  From  top  to  bottom:  density,  upward  mean 
velocity,  perpendicular  mean  energy,  parallel  mean  energy,  and  apex.  The  colours  are  associated  with  different  values  of  the  residence  time 
as  listed  in  the  toper  plots. 


be  simplified  by  introducing  the  following  approximation: 
Qi  ^  UiWi,  which  physically  means  that  w?  and  Wi  have 
the  same  scaling  versus  s  (Crew  et  al.,  1990;  Jasperse  and 
Grossbard,  2000).  Then,  Eq.  (23)  may  be  rewritten  as: 

(24) 

tD  oX  0*3 

where  w*  —  Wi  +  qi4>i  is  the  mean  total  energy.  When  in¬ 
tegrating  Eq.  (24)  over  the  box  limited  by  0  <  X'  <  X  and 
So  <  s'  <  s,  we  obtain  the  global  conservation  law  of  energy 


flux  densities  for  our  model: 

rs'^s 


r  '  Ni(X,s')w*(X,s')t^'dS'+ 

J  s'=so 

/  Ni(X',s)wJ(X',s)Ui(X\s)dX' 

t/x'=0 

Js'=so 

I  NiiX',  ioX(X'.  so)ui{X',  so)dX' 

Jx'=0 
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Ni(X\s')Wu(X\s')dX'ds\ 


(25) 


This  law  tells  us  how  much  energy  density  rate  goes  into  (first 
two  terms)  and  goes  out  of  (next  two  terms)  the  box,  limited 
by  0  <  X'  <  X  and  jq  <  <  “S*  This  is  balanced  by  the 

transferred  power  from  waves  to  particles  inside  the  same 
domain  (right-hand  side  term). 

Figure  5  displays,  for  the  previous  simulation,  the  varia¬ 
tions  of  the  outflow  rate  and  the  energy  density  rate,  as  fiinc- 
tions  of  s,  through  the  top  and  the  poleward  border  of  a  box 
limited  by  0  <  X  <  1 ,  <  •?'  <  -y-  As  shown  in  Fig.  5a,  the 

sum  of  the  two  contributions  is  still  equal  to  the  outflow  rate 
of  ions  injected  at  ^  =  sq,  and  therefore  Eq.  (23)  is  checked. 
As  shown  in  Fig.  5b,  the  conservation  of  the  energy  density 
flux  is  approximately  verified,  with  the  calculation  of  mo¬ 
ments  of  high  order  being  less  accurate  numerically. 


3.3  Sensitivity  of  the  results  to  the  boundary  conditions 


3.3.1  Influence  of  the  form  of  the  initial  distribution  fiinc- 
tion 


We  study  here  the  sensitivity  of  the  characteristics  of  the  ini¬ 
tial  ion  source  distribution  5o  on  the  steady-state  solution  of 
the  kinetic  equation.  Figure  6  shows  the  profiles  of  mo¬ 
ments,  along  the  poleward  heating  boundary  (X  =  1),  with 
the  parameters  wq  =  62  e V  and  a  =  1 ,  when  introducing  a 
parallel  temperature  in  the  source  distribution  Sq.  The  inter¬ 
esting  result  is  the  insensitivity  of  the  moment  profiles  to  the 
source  distribution  Sq.  Basically,  the  conjugate  action  of  the 
diffusion  and  mirror  forces  is  sufficiently  strong  so  that  the 
form  of  5*0  is  rapidly  rearranged  into  a  heated  distribution. 
For  ions  (not  shown),  the  effect  of  gravity  is  negligible 
and  the  change  on  moment  profile  is  therefore  insignificant. 

3.3.2  Influence  of  the  low-altitude  boundary 

By  considering  the  altitude  range  5-3  Re)  covered  by 
Interball-2  in  the  dayside  auroral  zone,  it  might  be  interesting 
to  analyze  the  sensitivity  of  the  results  to  the  lower  geocen¬ 
tric  distance  boundary,  sq,  of  the  simulation  domain,  where 
ionospheric  particles  are  injected.  Figure  7  shows  the  pro¬ 
files  of  O'*"  moments  along  the  poleward  heating  boundary 
(X=l)  when  considering  two  different  values  of  5*0 .  In  both 
cases,  the  heating  profile  is  the  same  as  used  above.  It  turns 
out  that  the  moments  are  more  slightly  sensitive  to  the  value 
of  So  than  the  characteristics  of  the  source  distribution  5o. 
However,  the  difference  in  upward  velocity  and  mean  energy 
components  is  small,  about  5%  and  ~  10%,  respectively. 
In  addition,  there  are  no  changes  in  the  asymptotic  evolution 
of  all  moment  profiles  versus  s. 


transversely  accelerated  ions  at  the  equatorward  edge  of  the 
dayside  cusp/cleft  (Andre  et  al.,  1988;  Whalen  et  al.,  1991). 
Here,  ions  may  be  preheated  before  entering  the  energization 
region,  which  results  in  significant  density  of  ions  at  high 
altitudes.  The  contribution  of  this  component  to  the  high- 
altitude  observations  usually  depends  on  the  residence  time 
to  and  the  strength  of  the  heating.  As  mentioned  in  Sect.  2.1, 
we  should  also  remind  the  reader  that  when  considering  a  ho¬ 
mogeneous  heating  region  in  latitude  and  a  constant  drift,  the 
ions  associated  with  the  population  flowing  from  the  topside 
ionosphere  and  the  one  from  the  equatorward  boundary  do 
not  mix.  Here,  for  the  values  of  considered,  we  found  that 
the  equatorward  population  leads  to  a  more  accurate  descrip¬ 
tion  of  the  low-latitude  side  of  the  heating  region  and  does 
not  drift  into  the  polar  cap  (X  >  1).  As  reported  in  earlier 
works  using  a  2-D  model  (Andre  et  al.,  1990),  the  results 
at  high  altitude  are  not  very  sensitive  to  the  temperature  of 
the  initial  distribution.  In  the  companion  paper  (Bouhram  et 
al.,  2003),  when  applying  the  results  to  a  few  specific  data 
sets,  the  contribution  of  the  equatorial  population  is  quanti¬ 
fied.  For  doing  so,  the  ion  density  may  be  inferred  from 
particle  measurements  at  the  satellite  altitude  Zsat-  Then, 
because  energization  processes  are  less  important  equator- 
ward  from  the  heating  region,  it  is  reasonable  to  assume  that 
the  density  is  inversely  proportional  to  the  cross-section  area 
A(5)  =  (s/so)^  of  a  magnetic  flux  tube  along  the  boundary, 
in  order  to  satisfy  the  particle  conservation  law. 

The  results  discussed  in  Sects.  3.3.1  and  3.3.2  point  out 
that  the  achieved  energies  of  ionospheric  ions  are  sufficiently 
important  at  higher  altitudes,  so  that  the  characteristics  of 
the  ion  source  distribution  have  been  forgotten,  as  reported 
in  previous  works  (Retterer  et  al.,  1987;  Crew  et  al.,  1990). 
In  the  same  way,  the  low-altitude  boundary  of  the  simula¬ 
tion  domain  has  a  weak  influence  on  the  results  at  higher 
altitudes.  This  justifies  that,  in  our  case  study,  an  accurate 
modelling  of  physical  processes  taking  place  at  low  altitude 
(<  2000  km)  is  not  necessary. 


4  Parametric  study 

4. 1  Effect  of  the  horizontal  drift 

In  this  section,  we  quantify  the  effects  of  the  parameter  tp, 
which  characterizes  the  horizontal  transport  in  the  poleward 
X  direction,  on  the  spatial  structure  of  ion  outflows.  Let  us 
remember  that  when  the  width  A  of  the  heating  region  is 
fixed,  to  is  proportional  to  the  poleward  convection  drift  Vp, 
Inside  the  heating  region,  we  also  analyze  the  differences  in 
terms  of  ion  energization  with  a  1-D  model,  where  the  con¬ 
vection  drift  is  turned  off  (Jd  oo). 


3.3.3  Influence  of  the  low-latitude  boundary 


4.1.1  Influence  inside  the  heating  region 


Another  possible  sensitive  parameter  is  the  distribution  of 
ions  along  the  equatorward  border  of  the  heating  region,  from 
which  ions  flow  into  the  region  under  the  effect  of  the  magne- 
tospheric  convection.  A  few  papers  reported  the  presence  of 


Figure  8  compares  H"^  and  O"*"  moment  profiles  along  the 
poleward  heating  boundary  (X  =  1)  for  different  values  of 
to  along  with  the  same  heating  profile  (iuq  =  62  eV,  a  =  1). 
Basically,  an  ion  drifts  out  of  the  heating  region  at  an  altitude 


1764 


M.  Bouhram  et  al.:  Modeling  transverse  heating  and  outflow  of  ionospheric  ions 


X 


Fig.  9.  Evolution  of  0“^  moments  versus  X  at  an  altitude  of  2  (a)  density,  (b)  upward  mean  velocity,  (c)  mean  energy,  (d)  fluence,  (e) 

perpendicular  mean  energy,  (f)  parallel  mean  energy,  (g)  apex,  (h)  energy  density  flux.  The  moments  associated  with  panels  (a),  (d)  and  (h) 
are  normalized  with  respect  to  the  values  at  the  low-altitude  boundaiy.  The  colours  are  associated  with  different  values  of  the  residence  time 
as  listed  in  panel  (a). 


Table  1.  Values  of  the  transferred  power  per  unit  area  and  flux  from  the  source  (at  1000  km)  from  waves  to  O"*"  and  ions,  as  a  function 

of  to 


to  (^) 

Uh^h)/ (•/H('So))(eV) 

200 

87 

160 

300 

139 

245 

500 

358 

315 

800 

622 

354 

+00 

940 

375 

proportional  to  its  resident  time  (<  to)  inside  the  heating  re¬ 
gion.  Therefore,  density  profiles  of  n,(5,  X  =  1)  decrease 
strongly  as  to  decreases.  In  the  same  way,  a  more  signifi¬ 
cant  fraction  of  ions  is  spread  out  of  the  heating  region 
(in  the  polar  cap)  for  smaller  f/).  For  time-of-flight  effects, 


upward  velocity  profiles  m/(5,  X  =  1)  increase  strongly  as  to 
decreases:  an  ion  reaching  a  fixed  distance  s  needs  a  larger 
upward  velocity  as  its  residence  time  decreases.  Because  of 
such  velocity  filter  effects,  ion  distribution  functions  have  a 
smaller  parallel  temperature  than  their  directed  parallel  en- 
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Fig.  10.  Evolution  of  H"*"  moments  versus  X  at  an  altitude  of  2  Re  >  The  format  is  the  same  as  in  Fig.  8. 


ergy  (miU^l),  Hence,  the  mean  parallel  energy  is  close 
to  (niiujll),  and  its  profile  u;,||(5,  X  =  1)  increases  strongly 
as  to  decreases.  On  the  other  hand,  perpendicular  tempera¬ 
ture  profiles  X  =  1)  are  insensitive  to  For  out¬ 

flowing  ions,  the  mirror  force  converts  the  perpendicular  en¬ 
ergy  into  parallel  energy  along  the  geomagnetic  field  lines. 
Therefore,  the  achieved  perpendicular  energies  are  only  due 
to  ICR  heating  and  do  not  depend  on  The  apex  angle  4'/, 
which  is  a  function  of  the  ratio  iu/x/u;/|j,  reaches  different 
asymptotic  values.  As  expected,  when  to  tends  to  infinity, 
tends  to  a  value  (43.3®)  as  given  by  1-D  models  (Crew  et 
al.,  1990): 


^  /w;/x\V2 

tan'I/i  =  ( — -)  = 

Viy,-  II / 


,iy/2  (6a  +  2)'/2 

.11/  ■  3 


Qualitatively,  when  to  is  larger  than  the  average  time  of  flight 
for  reaching  the  top  of  the  simulation  domain  (2  Re\  a  sig¬ 


nificant  part  of  the  ion  outflow  is  still  contained  inside  the 
heating  region.  Therefore,  all  moment  profiles  have  the  same 
shape  as  given  by  a  1-D  model,  on  which  the  ion  flux  tube  is 
completely  contained  inside  the  heating  region.  For  H"*"  ions, 
the  parameter  to  does  not  affect  the  moment  profiles.  Be¬ 
cause  the  H"*"  mass  differs  by  an  important  factor  (=16)  from 
the  0+  mass,  their  average  time  of  flight  through  the  simu¬ 
lation  domain  remains  always  smaller  than  typical  values  of 
to^ 

4.1.2  Influence  outside  the  heating  region 

Figures  9  and  10  compare  the  0+  and  H+  moment  profiles 
versus  X  at  .si  =  3  /?£  for  the  same  set  of  simulation  param¬ 
eters  as  in  Fig.  8.  We  note  that  the  perpendicular  temperature 
Wi±  and  apex  angle  profiles  fall  down  poleward  of  the 
heating  region  (X  >  1).  As  previously  mentioned,  an  ion 
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Fig.  11.  Evolution  of  0+  (left)  and  (right)  moments  along  the  poleward  heating  boundary.  From  top  to  bottom:  density,  upward  mean 
velocity,  perpendicular  mean  energy,  and  apex.  The  solid,  dashed  and  dotted  curves  are  associated  with  different  values  of  the  parameter  wq 
as  listed  in  the  toper  plots. 


drifting  out  of  the  heating  region  is  only  subject  to  the  mir¬ 
ror  force  and  gravity.  Because  of  the  absence  of  transverse 
energization,  heated  ion  distributions  are  rapidly  folded  into 
upwelling  ion  distributions  inside  a  thin  horizontal  region  as¬ 
sociated  with  the  sharp  gradient  along  X  of  Wij_  and 

Poleward  of  the  heating  region,  ions  are  filtered  owing  to 
their  time  of  flight  for  reaching  =  3  Re,  which  is  inversely 
proportional  to  their  mean  upward  drift  Uiis\,X).  We  will 
see  in  the  next  section  that  the  dispersion  Ui{s\ ,  X)  is  mainly 
governed  by  the  parameter  For  X  >  2,  the  apex  angle  'P/ 
of  increases  because  of  the  conjugate  action  of  poleward 
drift  and  gravity.  When  X  increases,  the  field-aligned  flow 
turns  downward  by  increasing  in  magnitude,  and  tends  to 
180°.  This  asymptotic  limit  is  not  observed  in  Fig.  9,  because 
the  simulation  domain  is  not  sufficiently  extended  in  X.  ince 
the  effect  of  gravity  is  negligible  for  'P/  does  not  increase 
for  X  >  2. 


The  effect  of  the  horizontal  drift  has  a  major  influence  on 
the  ionospheric  mass  loading  outside  the  heating  region  (po¬ 
lar  cap),  as  evidenced  in  density  ni(s\ ,  X),  fluence  Jiis\ ,  X), 
and  energy  density  flux  JiWi  (s\y  X)  profiles.  By  integrat¬ 
ing  Ji(s\,X)  and  JiWi{s\,  X)  over  X,  we  estimate  the  net 
outflow  rate  and  the  transferred  power  due  to  WPI  below 
s\  (Eqs.  22  and  25),  respectively.  For  smaller  values  of  to, 
the  energies  associated  with  the  maximum  fluence  decrease 
strongly  for  O"*".  In  the  same  way,  the  transferred  power  is 
as  small  as  to  is  small,  as  given  by  the  area  below  the  curves 
Ji  Wi  (5i ,  X).  Table  1  gives  the  transferred  power  from  waves 
to  and  H"’"  for  different  values  of  to-  We  have  added  the 
values  as  deduced  from  a  1-D  model  (to  — ►  -boo).  It  turns 
out  that  using  a  1-D  model  leads  to  a  strong  overestimation 
of  global  energy  density  transfers  to  heavy  ions,  such  as  0“^. 
On  the  other  hand,  the  effect  of  the  horizontal  drift  is  less  sig¬ 
nificant  for  H'*"  global  energy  density  transfers.  Once  again, 
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Fig.  12.  Evolution  of  0+  (left)  and  H"*"  (right)  apex  (top)  and  mean  perpendicular  energy  (bottom)  as  functions  of  wq  at  the  poleward  heating 
boundary  and  at  an  altitude  of  3  /?£.  Squares  and  diamonds  are  associated  with  a  spectral  index  a  =  1  and  or  =  2,  respectively.  The  solid 
curves  are  associated  with  the  analytic  formulas  from  the  1-D  theory. 


Fig.  13.  Evolution  of  the  dispersive  patterns  of  ion  upward  velocities  versus  X  at  two  altitudes:  2  Re  (left)  and  3  Re  (right).  The  curves  are 
displayed  from  tp  =  200  s  up  to  tp  =  800  s  with  a  step  of  100  s. 


this  difference  is  due  to  their  lower  mass  and  leads  to  average 
times  of  flight  smaller  than  tp, 

4.2  Influence  of  the  heating  profile 

Figure  1 1  compares  the  0"*“  and  }l^  moment  profiles  versus 
X  at  =  4/?£  for  different  heating  profiles  with  a  =  1 
fixed,  and  tp  =  500  s.  For  a  fixed  value  of  tp,  the  velocity 
dispersion  in  X  remains  insensitive  to  the  heating  parame¬ 
ters.  On  the  other  hand,  Figs.  9  and  10  have  pointed  out  that 
the  dispersion  Ui(s,  X)  exhibits  different  shapes  when  con¬ 
sidering  different  values  of  tp.  Therefore,  we  may  conclude 
that  the  simulated  velocity  dispersion  versus  X  is  a  function 
of  tp.  This  result  is  important  when  studying  an  upfl owing 
ion  event,  because  it  suggests  that  the  convection  drift  can 
be  determined  from  the  observed  velocity  filter  effect.  Let  us 
introduce  the  parameters  Wi±p  and  '1^/^,  corresponding  to  the 


mean  perpendicular  energy  and  the  apex,  respectively,  at  the 
poleward  heating  boundary  (PHB)  at  a  distance  s  for  a  type 
of  ion  i.  As  shown  in  Fig.  11,  the  value  of  WiXp  increases 
strongly  versus  wq-  Conversely,  the  value  of  ^ip  is  approxi¬ 
mately  constant  for  H"*"  and  slightly  dispersed  for  O'*"  (~5°). 
In  the  following,  we  will  see  that  Wi±p  and  ^ip  may  be  used 
as  a  probe  of  the  heating  parameters  (u;o,  of). 

Figure  12  displays  the  variation  of  ^tp  and  Wij_p  at  s  = 
4  Re  as  functions  of  wq  and  a.  Solid  curves  correspond  to 
the  variations  of  Wi±p  and  ^ip  in  the  1-D  case.  In  this  ap¬ 
proximation,  '1'/^  is  given  by  Eq.  (26)  and  depends  on  a  only. 
For  Wi±p,  there  is  no  analytic  formula.  However,  on  the  basis 
previous  studies  (Crew  et  al.,  1990;  Barghouthi,  1997),  Wi±p 
may  be  approximated  as  the  perpendicular  energy  Wi±^  of  a 
single  mean  ion  (Eq.  20).  As  evidenced  in  Fig.  12,  the  val¬ 
ues  of  WiXp  remain  close  to  the  1-D  limit  given  by  Eq.  (20). 
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Fig.  14.  Patterns  in  a  plane  luo  —  « 
of  the  apex  and  the  mean  perpendicu¬ 
lar  energy  of  ions  at  the  poleward 
edge  of  the  heating  region  at  two  alti¬ 
tudes:  1.5 /?£  (left)  and  3  Re  (right). 
The  oblique  curves  correspond  to  con¬ 
tours  of  the  mean  perpendicular  energy 
(in  eV)  while  the  curves  approximately 
horizontal  correspond  to  contours  of  the 
apex  with  a  constant  step  of  2°.  Crosses 
correspond  to  values  inferred  from  nu¬ 
merical  simulations. 


Slight  differences  appear  on  ^ip  profiles  for  where  val¬ 
ues  remain  lower  than  the  1-D  limit,  because  of  limited  res¬ 
idence  times  (see  Sect.  4.1.1).  However,  when  a  increases, 
^ip  changes  little  and  tends  to  the  1-D  limit. 

4.3  Parametric  survey 

In  the  last  section,  we  have  studied  qualitatively  how  the  sim¬ 
ulation  parameters  (tp,  wo,  a)  control  the  moment  evolution 
at  higher  altitudes.  Furthermore,  we  have  established  that  in 
the  PHB  the  mean  perpendicular  energy  wup  and  the  apex 
^ip  associated  with  high-altitude  heated  conics  are  mainly 
controlled  by  the  heating  parameters,  wq  and  a.  Conversely, 
outside  the  heating  region,  the  slope  dLog(M/)/dLog(X)  char¬ 
acterizing  the  velocity  filter  effect  is  insensitive  to  the  heating 
parameters  and  depends  on  tp. 

In  this  section,  by  gathering  the  results  of  a  large  set 
of  simulation  runs,  we  construct  a  unique  relationship  be¬ 


tween  the  triplet  [wi±p,  ^ip,  dLog(ui)/dLog(X)],  associated 
with  parameters  which  are  experimentally  available,  and  the 
triplet  (fp.  Wo,  a)  associated  with  the  horizontal  drift  and  the 
transverse  heating. 


4.3.1  Evolution  of  the  velocity  filter 

Figure  13  shows  the  evolution  of  the  ion  dispersion 
Log[Mi(X)]  for  different  values  of  tp  at  two  geocentric  dis¬ 
tances.  Here,  the  slope  dLog(M|)/dLog(X)  turns  out  to  be 
flatter  at  higher  altitudes.  Let  us  consider  the  flow  line  u\ 
crossing  the  point  (^i,  Xi).  From  Eq.  (21),  at  a  geocentric 
distance  52  =  Jfi  +  A^,  the  line  u\  is  translated  along  X  by 
a  distance  AX  =  x  Aj/mi,  inversely  proportional  to 
Ml.  Therefore,  the  slope  Log[M,(X)]  decreases  flatly  versus 
X  between  si  and  52. 
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Fig.  15.  Patterns  in  a  plane  i^o  —  of 
of  the  apex  and  the  mean  perpendicular 
energy  of  ions.  The  format  is  the 
same  as  in  Fig.  14. 


4.3.2  Relationship  between  heating  and  conic  parameters 


From  a  large  set  of  simulations,  we  may  construct,  at  a  fixed 
altitude,  contour  plots  of  Wi±p  and  as  functions  of  the 
ICR  parameters  (luo,  a).  Figures  14  and  15  provide  contour 
plots  of  Wi±p  and  ^ip  for  and  H"^,  respectively,  by  con¬ 
sidering  different  altitudes  and  values  of  tp-  These  contours 
have  been  constructed  from  1-D  fits,  similar  to  Fig.  12.  In 
the  limit  of  higher  u;o  or  a,  the  effect  of  the  horizontal  drift 
becomes  negligible  and  the  contours  follows  those  given  by 
the  1-D  theory  (Eqs.  20  and  26).  In  the  domain  (wq  <  1  eV, 
a  <  1),  contours  are  not  shown  because  such  heating  param¬ 
eters  lead  to  ion  distributions  associated  with  insignificant 
energies  (lu/xp  <  10  eV)  for  being  observed  at  the  PHB.  At 
a  fixed  altitude,  the  Wi±p  curves  have  flatter  slope  for  0“^ 
than  By  differentiating  Eq.  (20)  and  retaining  only  dom¬ 


inant  terms,  the  slope  can  be  expressed  as: 

r  1  ^ _ \ _  (27) 

LaLog(u)o)J„,,j.  Log[M/^^(s/so)^] 

Therefore,  the  slope  of  lu/xp  contour  decreases  with  respect 
to  the  ion  mass  and  the  altitude.  Furthermore,  there  is  a 
single  intersection  between  a  lu/xp  contour  and  a  'Pip  con¬ 
tour.  Therefore,  a  couple  {wq,  a)  leads  at  a  fixed  altitude  to  a 
unique  ion  conic  form  at  the  PHB  and  is  characterized  by  the 
couple  of  moments  {wi±p,  ^ip). 

5  Summary  and  conclusion 

In  this  paper,  we  have  studied  the  spatial  structure,  from 
1000  km  up  to  20  000  km  in  altitude,  of  dayside  ionospheric 
ion  outflow  as  a  function  of  ion  heating  parameters  and 
the  convection  drift  using  a  two-dimensional  (2-D),  Monte 
Carlo,  trajectory-based  code.  Since  no  analytic  theory  is 
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available  in  2-D,  the  results  of  the  model  have  been  tested 
successfiilly  on  the  basis  of  global  conservation  laws. 

In  contrast  to  previous  work  based  on  1-D  simulation,  the 
introduction  of  a  poleward  convection  drift  leads  to  a  lim¬ 
ited  residence  time  Id  of  ions  when  being  heated  in  the  day- 
side  cusp/cleft.  Residence  time  effects  turn  out  to  be  impor¬ 
tant  by  modifying  the  density  patterns  of  O'*"  outflows,  while 
residence  time  effects  on  density  patterns  are  small  for  H"*" 
outflows.  These  effects  also  control  the  dispersive  energy 
patterns  of  outflows  poleward  of  the  heating  region  (i.e.  in 
the  polar  cap).  Conversely,  inside  the  heating  region,  lim¬ 
ited  residence  times  lead  to  moderate  and  weak  effects  in  the 
angle  and  energy  characteristics  associated  with  O'*"  and  H'*' 
conics,  respectively.  Such  parameters  are  to  the  first  order 
controlled  by  the  strength  of  ion  transverse  heating,  as  char¬ 
acterized  in  the  model  by  the  couple  of  parameters  {wq,  a). 
In  terms  of  energy  density  flux  exchanges,  residence  time 
effects  are  important  for  heavy  populations,  such  as  O"*".  Be¬ 
cause  the  ion  density  is  overestimated  using  a  1-D  model,  the 
introduction  of  a  limited  residence  time  tp  leads  to  smaller 
transfers  of  power  between  the  waves  and  the  O'*"  ions  inside 
the  heating  region. 

The  main  conclusion  of  the  paper  is  that  the  transport  pat¬ 
terns  of  ion  outflows  when  associated  with  a  reasonable  num¬ 
ber  of  free  parameters  u;o,  a),  lead  to  unique  features  in 

space.  In  a  companion  paper  (Bouhram  et  al.,  2003),  pro¬ 
viding  that  our  assumptions  are  satisfied,  we  reverse  the  pro¬ 
cess  by  determining  the  strength  of  ion  heating  (i^o,  or)  for  a 
large  number  of  data  sets  in  the  dayside  cusp/cleft  from  the 
Interball-2  satellite. 
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